class ASTroCAT_Chemostat extends Environment {
  float seed = 0.01;
  ASTroCAT eco;
  
  ASTroCAT_Chemostat(int NP) {this(NP,1);}
  
  ASTroCAT_Chemostat(int NP, int NPC) {
    eco = new ASTroCAT_Clonal(NP,NPC);
    linkToEcosystem(eco);
    for (int j=0; j<NP; j++) {
      getVar("P-" + j).setInitial(seed/NP);
      getVar("Z-" + j).setInitial(seed/NP);
    }
    
    // chemostat params and fluxes
    addParam("external N concentration", "N0", "µM", 30, 0, 100);
    addParam("relative flowthrough rate", "q", "1/d", 0.025, 0, 0.2);
    addParam("flowthrough intermittency", "n_imt", "", 0, 0, 20).quantize(1);
    addParam("flowthrough period", "T", "d", 10, 2, 365);
    addFlux("N supply", "supply", null, eco.N, 1);
    for (int i=0; i<eco.vars.length; i++) {
      String theVar = eco.vars[i].shortname;
      if (eco.vars[i].advects && eco.vars[i] != eco.N) {
        addFlux(eco.vars[i].shortname + " flowthrough", "flow_" + eco.vars[i].shortname, eco.vars[i], null, 1);
      }
    }
    
    // light params
    addParam("depth", "H", "m", 20, 1, 100).logScale = true;
    addParam("surface PAR", "PAR0", "W/m^2", 150, 0, 300);
    addParam("PAR seasonal cycle on/off", "do_PARcycle", "", 0, 0, 1).quantize(1);
    
    // sinking
    addParam("SD sinking rate", "wsinkSD", "m/day", 10, 0, 100);
    addParam("LD sinking rate", "wsinkLD", "m/day", 100, 0, 100);
    addFlux("SD sinking", "sink_SD", eco.SD, null, 1);
    addFlux("LD sinking", "sink_LD", eco.LD, null, 1);

    startArchive(1);
  }


  void calcFluxes() {
    // calc light
    getVar("PAR").current[0] = getParamVal("PAR0") * exp(-0.5 * getParamVal("H") * (getParamVal("attSW") + getParamVal("attP") * getVar("Ptot").current[0]));
    if (getParamVal("do_PARcycle")==1) getVar("PAR").current[0] *= 0.5 * (1 - cos(2*PI*getVar("t").current[0]/365.));
    // calc chemostat flowthrough
    getFlux("supply").current[0] = getParamVal("q") * (getParamVal("N0") - eco.N.current[0]);
    if (getParamVal("n_imt") > 0) {
      getFlux("supply").current[0] *= pow(sin(2*PI*getVar("t").current[0]/getParamVal("T")), 2*getParamVal("n_imt"));
    }
    for (int i=0; i<eco.vars.length; i++) {
      Flux F = getFlux("flow_" + eco.vars[i].shortname);
      if (F != null) {
        F.current[0] = getParamVal("q") * eco.vars[i].current[0];
      }
    }
    // sinking
    getFlux("sink_SD").current[0] = getParamVal("wsinkSD") / getParamVal("H") * eco.SD.current[0];
    getFlux("sink_LD").current[0] = getParamVal("wsinkLD") / getParamVal("H") * eco.LD.current[0];
    eco.calcFluxes();
  }
  
  
  void saveArchive() {
//    archive.saveFluxes = false;
    ((ASTroCAT)eco).saveSpecialInfo(archive);
    super.saveArchive();
  }
  
}



// ---------------------------------------------------------------------------------------------------



class ASTroCAT_Upwelling extends Slice {

  ASTroCAT_Upwelling(int M, int K, float Lshelf, float Ltot, int NP, float Ntot) {
    name = "2D upwelling conveyor";
    shortname = "upwelling";
    
    float Hwall = 30, Hshelf = 130, stretching = 8;
    
    dx = new float[M];
    for (int m=0; m<M; m++) dx[m] = Ltot/M;
    dz = new float[K];
    float dz0 = 2 * Hshelf / K / (stretching+1);
    for (int k=0; k<K; k++) dz[k] = dz0 * (1 + (stretching-1)/(K-1) * (K-1-k));
    createDomain(dx,dz);
    float[] depth = new float[M];
    float slope = (Hshelf - Hwall) / Lshelf;
    for (int m=0; m<M; m++) {
      depth[m] = slope * (x_rho[m] - Ltot) - Hwall;
    }
    setDepths(depth);
    
    linkToEcosystem(new ASTroCAT(NP));
    eco.getParam("mu0").hidden = true;
    eco.getParam("I0").hidden = true;
    eco.getParam("K0").hidden = true;
    eco.getParam("xpreyopt").hidden = true;
    eco.getParam("epsilon").hidden = true;
    eco.getParam("maxPref").hidden = true;
    
    
    float seed = 0.1;
    for (int j=0; j<NP; j++) {
      getVar("P-" + j).setInitial(seed/NP);
      getVar("Z-" + j).setInitial(seed/NP);
    }
    getVar("N").setInitial(max(Ntot-2*seed, 0));
    getVar("SD").setInitial(0);
    getVar("LD").setInitial(0);
    
    internalTimestep = 0.001;

    Statevar mass = addVar("mass","mass","",Ncells);
    getVar("mass").setInitial(1);
    mass.advects = true;
    mass.mixes = true; 
    Statevar age = addVar("age","age","d",Ncells);
    age.advects = true;
    age.mixes = true;
    addFlux("aging","aging",null,age,Ncells);
    
    leftBC = new BoundaryCondition(this);
    leftBC.addVar("mass",1);    
    leftBC.addVar("age",0);
    
    allocateStorageForPhysics();

    addParam("surface PAR","PAR0","W/m2", 150, 0, 300);
    addParam("log vertical mixing","logkappa","m2/s", -4.5, -6, -4);
    float lossFrac = 1;
    addParam("SD sinking rate", "wsinkSD", "m/day", 1, 0, 100);
    addSinking("SD", getParamVal("wsinkSD"), lossFrac);
    addParam("LD sinking rate", "wsinkLD", "m/day", 100, 0, 100);
    addSinking("LD", getParamVal("wsinkLD"), lossFrac);
    
    addParam("mean surface velocity", "u00", "m/s", -0.05, -0.3, 0.3);
    addParam("Ekman layer thickness", "hsurf", "m", 20, 0, 40);
    addParam("surface velocity variance", "du0", "m/s", 0, 0, 0.3);
    addParam("surface velocity period", "Tu", "d", 10, 2, 365);
    startArchive(1);
  }
  
  float attTot(int k, int m) {
    float att = getParamVal("attSW");
    att += getParamVal("attP") * getVar("Ptot").current[km2ind(k,m)];
    return att;
  }

  void updatePrivateParams() {
    float kappa0 = pow(10,getParamVal("logkappa")) * 86400;
    for (int m=0; m<M; m++) for (int k=0; k<K; k++) kappa[k][m] = kappa0;
    sinkingRate[getVarIndex("SD")] = getParamVal("wsinkSD");
    sinkingRate[getVarIndex("LD")] = getParamVal("wsinkLD");
    super.updatePrivateParams();
  }
  
  void calcFluxes() {
    for (int m=0; m<M; m++) {
      for (int k=0; k<K; k++) {
        getFlux("aging").current[km2ind(k,m)] = 1;
      }
    }
    super.calcFluxes();
  }
  
  void calcVelocity() {
    float u0 = 86400 * getParamVal("u00") + 86400 * getParamVal("du0") * sin(TWO_PI * getVar("t").current[0] / getParamVal("Tu"));
    float hsurf = getParamVal("hsurf");
    for (int m=0; m<M+1; m++) for (int k=0; k<K; k++) u[k][m] = 0;
    // apply velocity in the surface layer
    for (int m=0; m<M; m++) {
      for (int k=kbotu[m]; k<K; k++) {
        u[k][m] = u0 * exp(z_rho[k]/hsurf);
      }
    }
    // apply a compensating barotropic velocity
    for (int m=0; m<M+1; m++) {
      float uAtot = 0, Atot = 0;
      for (int k=kbotu[m]; k<K; k++) {
        uAtot += u[k][m] * Ayz[k][m];
        Atot += Ayz[k][m];
      }
      float umean = uAtot / Atot;
      for (int k=kbotu[m]; k<K; k++) u[k][m] -= umean;
    }
    super.calcVelocity();
  }
  
  void saveArchive() {
    ((ASTroCAT)eco).saveSpecialInfo(archive);
    super.saveArchive();
  }

  
}




// -------------------------------------------------------------------------------------------------------------




class ASTroCAT_WaterColumn extends Slice {

  ASTroCAT_WaterColumn(int K, float Htot, int NP, float Ntot) {
    name = "1D water column";
    shortname = "watercol";
    createDomain(1, K, 1, Htot);
    linkToEcosystem(new ASTroCAT(NP));
    
    float seed = 0.01;
    for (int j=0; j<NP; j++) {
      getVar("P-" + j).setInitial(seed/NP);
      getVar("Z-" + j).setInitial(seed/NP);
    }
    getVar("N").setInitial(max(Ntot-2*seed, 0));
    getVar("SD").setInitial(0);
    getVar("LD").setInitial(0);

    internalTimestep = 0.001;

    Statevar mass = addVar("mass","mass","",Ncells);
    getVar("mass").setInitial(1);
    mass.advects = true;
    mass.mixes = true; 
    
    allocateStorageForPhysics();

    addParam("surface PAR","PAR0","W/m2", 150, 0, 300);
    addParam("vertical mixing","logkappa","m2/s", -5, -6, -4);
    
    addParam("SD sinking rate", "wsinkSD", "m/day", 1, 0, 100);
    addSinking("SD", getParamVal("wsinkSD"), 0);
    addParam("LD sinking rate", "wsinkLD", "m/day", 100, 0, 100);
    addSinking("LD", getParamVal("wsinkLD"), 0);

    startArchive(1);
  }
  
  float attTot(int k, int m) {
    float att = getParamVal("attSW");
    att += getParamVal("attP") * getVar("Ptot").current[km2ind(k,m)];
    return att;
  }

  void updatePrivateParams() {
    float kappa0 = pow(10,getParamVal("logkappa")) * 86400;
    for (int m=0; m<M; m++) for (int k=0; k<K; k++) kappa[k][m] = kappa0;
    sinkingRate[getVarIndex("SD")] = getParamVal("wsinkSD");
    sinkingRate[getVarIndex("LD")] = getParamVal("wsinkLD");
    super.updatePrivateParams();
  }
  
  void saveArchive() {
    ((ASTroCAT)eco).saveSpecialInfo(archive);
    super.saveArchive();
  }


}

